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Abstract. We report on a computational approach based on the self-consistent solution of the steady-state 
Boltzmann transport equation coupled with the Poisson equation for the study of inhomogeneous transport in deep 
submicron semiconductor structures. The nonlinear, coupled Poisson-Boltzmann system is solved numerically us- 
ing finite difference and relaxation methods. We demonstrate our method by calculating the high-temperature 
transport characteristics of an inhomogeneously doped submicron GaAs structure where the large and inho- 
mogeneous built-in fields produce an interesting fine structure in the high-energy tail of the electron velocity 
distribution, which in general is very far from a drifted-Maxwellian picture. 
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1. Introduction 

The carrier dynamics in submicron structures is far 
from thermal equilibrium due to strong and rapidly 
varying external and built-in electric fields. Hot elec- 
tron and ballistic effects dominate the transport char- 
acteristics and the electron velocity distribution func- 
tion in such systems is far from a drifted-Maxwellian 
description. In order to fully take into account the 
nonequilibrium nature of the transport, a full solu- 
tion of the semiclassical Boltzmann transport equation 
(BTE) is required. Although the Monte Carlo method 
has been very popular for the solution of the BTE in 
semiconductor device simulation 0, several works 
[7] have recently solved the BTE by direct methods, 
thus allowing noise-free spatial and temporal resolu- 
tion of the electron distribution function, which in the 
Monte Carlo method may be difficult to obtain due to 
the statistical nature of the approach. In this paper, 
we present a straight-forward approach to calculate the 
electron distribution function, f(x,v), for submicron 
inhomogeneous semiconductor structures by solving 
the steady-state BTE self-consistently with the Pois- 
son equation. We solve the strictly two-dimensional 



(2D) BTE (one dimension corresponding to position 
and one to velocity) and treat scattering within the 
relaxation time approximation (RTA) where each indi- 
vidual scattering mechanism is represented by a char- 
acteristic scattering rate that can be derived from 
quantum mechanical scattering theory. We demon- 
strate our approach for submicron, inhomogeneously 
doped structures and discuss the general nonequilib- 
rium transport characteristics. 

2. Basic equations 

The Boltzmann equation describes the dynamics of 
the semiclassical distribution function, /(r, v, t), un- 
der the influence of electric and magnetic fields, as 
well as different scattering processes. In the absence 
of a magnetic field, the 2D phase-space, steady-state 
BTE in the RTA is written according to: 



eE(x) df(x, v) ^ df(x, v) 



dv 



dx 



f(x,v) - Ile{x,v) 
r{e) 

(1) 

where m* is the electron effective mass in the parabolic 
band approximation, and fhE{x,v) is a local equilib- 
rium distribution function appropriate to a local den- 
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sity, applied field and equilibrium lattice temperature, 
To, to which the distribution function f(x, v) relaxes at 
a relaxation rate r(e)" 1 . As the local equilibrium func- 
tion, we choose in the following a Maxwell-Boltzmann 
(MB) distribution at To, normalized to the local den- 
sity n(x) 



Ile{x,v) = n(x) 



-l 1/2 



27rfcT f 



(2) 



The inhomogeneous electric field, E(x), in the BTE, 
originating from the spatially dependent electron and 
doping densities, n(x) and Nu(x), is obtained from 
the Poisson equation 



d 2 (t> _ _dE 
dx 2 dx 



e Nn{x)-n{xl = 



where e is the static dielectric constant. Since the elec- 
tron density is related to the distribution function by 



n(x) = I f(x, v)dv 



(4) 



the Poisson and Boltzmann equations constitute a cou- 
pled, nonlinear set of equations, and thus, Eqs. Ill II) 
have to be solved self-consistently. 

3. Numerical procedure 

The numerical procedure consists, in short, of initializ- 
ing the system parameters, discretizing Eqs. (Ill 1 fi on a 
2D grid in phase-space, performing the self-consistent 
Poisson-Boltzmann loop and, upon convergence, cal- 
culate and output the electron distribution function, 
electric field and the desired moments of the BTE. In 
the calculations, after initialization, we rescale the sys- 
tem parameters and the equations according to 



c' = x/Ld, v' 



vt/Li 



(5) 



where L D = ^ee kBT /e 2 N is the Debye length, 
N = ma,x[Ni>(x)] and r is a characteristic scattering 
time. The choice of grid size and resolution depends 
to a large extent on the system parameters and the 
electrostatics present in the device. In order to repro- 
duce details due to strong and rapidly varying electric 
fields, we choose the spatial grid step size to be smaller 
than the Debye length, Ljj, defined above. In velocity 
space, on the other hand, the discrete grid step size 
needs to be small enough to resolve fine structure in 
the distribution function, as well as give accurate re- 
sults for the moments of the BTE. In addition, the 



grid needs to be large enough, in velocity, in order to 
capture the full information in the high-energy tail of 
the distribution function, and in position, in order to 
damp out the effects of the contact boundaries. 

The Poisson and Boltzmann equations are solved by 
finite difference and iterative relaxation methods 8 . 
For the Poisson equation J3J, we use forward and back- 
ward Euler differences according to 



L X L X 4>j 



4>j+\ - ^4>j + • 
(5x) 2 



-Pj 



(6) 



where L+<fi(x) — (4>j+i — <t>j)/Sx and L~<f>(x) = (<fij — 
<j>j—i)/8x denote forward and backward Euler steps, 
respectively. The resulting matrix equation is solved 
iteratively using successive overrelaxation (SOR) [Sj. 

For the solution of the BTE, we adopt an upwind 
finite difference scheme 3 which amounts to the fol- 
lowing discretization of the partial derivatives in Eq. 



0/ 
dv 
Of 
dx 



= L+Mf(x,v) E(x) > [E{x) < 0] (7) 



L+Mf(x,v) v < [v > 0] 



(8) 



As for the Poisson equation, we use SOR to solve the 
matrix equation resulting from the discretization of 
Eq. 0. 

For the boundary conditions of the Poisson- 
Boltzmann equations we adopt the following: For the 
potential, the values at the system boundaries, de- 
noted (l)eft and (r)ight are fixed to 4>{x{) — Uq and 
4>{x r ) = 0, respectively, corresponding to an externally 
applied voltage Uq. The electron density is allowed to 
fluctuate freely around the boundaries, subject to the 
condition of global charge neutrality, which is enforced 
between each successive iteration in the self-consistent 
Poisson-Boltzmann loop. We choose the size of the 
highly-doped contacts to be large enough such that 
the electron density and the electric field deep inside 
the contacts is constant. 

For the electron distribution function four bound- 
ary conditions can be defined in the 2D phase-space. 
At the velocity cut-off in phase-space, we choose 
f(x,v rnax ) = f(x,-v max ) = fLE{x,v) which is rea- 
sonable since we assume v max > 30fcsTo in the cal- 
culations. At such high velocities, the electron pop- 
ulation is negligible and of the same order as the lo- 
cal equilibrium distribution /le(x,v). At the contact 
boundaries, we assume that the electric field is low 
and constant (as verified in the calculations), and thus, 
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the homogeneous solution to the BTE in the linear re- 
sponse regime of transport applies. Hence, 

f{xi,v) = f LE (xi,v)[l - vE{xi)T(e)/k B T ], (9) 

where i — l,r. The iterative Poisson-Boltzmann loop 
consists of an updating procedure for the electric field, 
electron distribution function and electron density us- 
ing Eqs. (0nEJl> until convergence. The convergence 
criterion is determined and checked in terms of the 
evolution of the Li norm of the potential and density 
variations between subsequent iterations. Typically, 
the results are converged when the Li norms for the 
potential and density are on the order of 10~ 3 of the 
original values. Between subsequent iterations, we em- 
ploy linear mixing in the electron density, according to 

n\x) = (1 - a)n old {x) + an new (x) , (10) 

where n old (x) is the input density to the Poisson 
solver, n new (x) is the new density obtained from the 
solution of the BTE using the new electric field ob- 
tained from the Poisson solver, and n'(x) is the final 
density that is used as an input to the next iteration 
in the Poisson-Boltzmann loop. The convergence and 
stability of the self-consistent loop are strongly depen- 
dent on the system parameters and the nonequilib- 
rium nature of the electronic system. If the system 
is strongly out of equilibrium, displaying large vari- 
ations and strengths of the electric field, the mixing- 
parameters a may have to be chosen as small as a few 
percent, thus affecting the overall runtime. Further- 
more, for highly doped structures, the convergence is 
slower, partly due to the required small grid size in 
position due to the small Debye length, but also due 
to the slow convergence in the SOR procedure in the 
BTE, where the stability of the numerical scheme is 
given in terms of a Courant-Friedrich-Levy type con- 
dition 8 . Still, the computational demands for the 
calculations reported in this paper are modest. 

4. Numerical results 

In the following, we demonstrate our numerical ap- 
proach with calculations of the transport character- 
istics of a model GaAs n + ~ n~ — n + — n~ — n + 
structure with the doping densities n + = 10 23 m~ 3 and 
n~=10 19 m~ 3 . In order to highlight the effects of in- 
homogeneities and scattering while keeping the nature 
of the scattering structureless, we use a constant scat- 
tering time t — 2.5T0" 13 s, which corresponds to real- 
istic mobilities of GaAs at room temperature for which 



the calculations have been performed. The central 
n~ — n + — n~ region has the dimensions 200/200/200 
nm, whereas the contacts are 1 /im long. 



n + , . n\ . _ n + , , . vc , n + 
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Figure 1: Electron potential energy and electric field in the 
central region of the n + — n~ — n + — n~ — n + system with 
dimensions described in the text. The corresponding current 
density is j(x) = —5.4 • 10 4 A/cm 2 . 

In Fig. ^we show the electric field and potential en- 
ergy around the central region of the system described 
above, subject to an applied bias voltage Vb = —0.5 V. 
Due to the charge imbalance, electrons diffuse towards 
the lightly doped regions, where potential barriers are 
formed and, correspondingly, a large and inhomoge- 
neous electric field on the order of 10 kV/cm is formed, 
even in the absence of an external applied voltage. As 
a finite voltage is applied to the device, the majority 
of the potential drop occurs over the submicron cen- 
tral region, giving rise to a strongly inhomogeneous 
field distribution, in contrast to the n + contact re- 
gions, where the field in comparison is very low and 
constant. 

The electron velocity distribution in the central re- 
gion of the structure is shown in Fig. Ufa), for five 
specific spatial points as depicted in Fig. ^ Figure 
Gib) shows a contour plot of the full spatial depen- 
dence of f(x,v) in that region. It is clear that the 
inhomogeneous electric field gives rise to a strong spa- 
tial dependence of the velocity distribution function 
along the direction of transport, and that the distri- 
bution function in the central region is very far from 
thermal equilibrium. 

In the outermost highly doped n + regions, where 
the field is low and constant the distribution is sim- 
ply a shifted Maxwellian. In the lightly doped n~ 
regions on the other hand, the velocity distribution is 
highly asymmetric and develops a narrow peak that 
rapidly shifts toward higher velocity along the direc- 
tion of transport. This peak contains quasi-ballistic 
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electrons which are accelerated by the strong electric 
field in the central region, and thus, have a consider- 
ably larger average velocity compared to the electrons 
in the contacts. Close to the potential barrier the dis- 
tribution function is suppressed at low velocities due to 
the skimming of the distribution of incoming electrons, 
as well as the restriction of drain induced electron flow 
with v < due to the potential barrier. 

However, the low-velocity contribution to the dis- 
tribution function gradually increases away from the 
barrier, as thermionically injected electrons gradu- 
ally are thermalized and the lower effective barrier 
height allows electrons from the n + regions to pen- 
etrate the lightly-doped region. Thus, the total dis- 
tribution function consists of a quasi-ballistic, high- 
velocity and a diffusive, low-velocity contributions, 
which gives the total distribution function a highly 
non-Maxwellian broad and asymmetric shape. Fur- 
thermore, the presence of the two barriers creates an 
additional quasi-ballistic structure in the high-energy 
tail of the distribution function in the second n~ re- 
gion, as electrons that traverse the intermediate n + re- 
gion ballistically get an additional acceleration toward 
higher velocities by the electric field in the second rT 
region, thus creating two high- velocity electron beams. 
These features emphasize the highly nonequilibrium 
nature of the electron transport in these type of sys- 
tems and demonstrate that our method is capable of 
taking them fully into account. 
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Figure 2: (Color online) Normalized electron velocity dis- 
tribution function f(x,v)v t h/n(x), where v t ^ is the ther- 
mal velocity, at (a) five different points in space correspond- 
ing to the points depicted in Fig. (b) Contour plot of 
f(x,v)v th /n(x). 



5. Conclusions 

We have presented a numerical method for the solution 
of the steady-state, coupled Poisson-Boltzmann equa- 
tions for the study of inhomogeneous, submicron semi- 
conductor structures and demonstrated our approach 
on a submicron GaAs structure with strong built-in 
electric fields. We have shown that our method is ca- 
pable of taking into account the strong nonequilibrium 
transport properties that arise in such systems due to 
the presence of very large and inhomogeneous electric 
fields, and that interesting structure is present in the 
high-energy tail of the distribution function, caused by 
quasi-ballistic electrons. 
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